Phys. Rev. E ?, ? (1998), accepted for publication. 



1 



Logarithmic corrections of the avalanche distributions of sandpile models at the upper 

critical dimension 

S. LlibeckQ 

Theoretische Tieftemperaturphysik, Gerhard-Mercator-Universitdt Duisburg, 
Lotharstr. 1, 47048 Duisburg, Germany 
(Received 20 May 1998) 

We study numerically the dynamical properties of the BTW model on a square lattice for various 
dimensions. The aim of this investigation is to determine the value of the upper critical dimension 
where the avalanche distributions are characterized by the mean-field exponents. Our results are 
consistent with the assumption that the scaling behavior of the four-dimensional BTW model is 
characterized by the mean-field exponents with additional logarithmic corrections. We benefit in 
our analysis from the exact solution of the directed BTW model at the upper critical dimension 
which allows to derive how logarithmic corrections affect the scaling behavior at the upper critical 
dimension. Similar logarithmic corrections forms fit the numerical data for the four- dimensional 
BTW model, strongly suggesting that the value of the upper critical dimension is four. 
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I. INTRODUCTION 

The concept of self-organized criticality introduced by 

[ Bak, Tang, and Wiesenfeld allows to describe scale invari- 
ance in driven systems 0. Sandpile models and espe- 

[ cially the Bak- Tang- Wiesenfeld (BTW) sandpile model 
are known as the paradigm of self-organized criticality. 
The steady state dynamics of the system is character- 
ized by the probability distributions for the occurrence 
of relaxation clusters of a certain size, area, duration, 

I etc. Despite numerous theoretical efforts the val- 

ues of the exponents of the probability distribution char- 

[ acterizing the critical behavior of the system were de- 
termined only numerically for D = 2 and D — 3 [^,0. 
These investigations are based on an accurate finite-size 
scaling analysis and were confirmed in a recently pub- 
lished work In higher dimensions the scaling behav- 
ior of the BTW model is still controversial. Especially, 
the value of the upper critical dimension where the 
mean-field solution describes the scaling behavior of the 
system is not known exactly. Whereas renormalization 
group approaches predicted D„ = 4 |p|-pd| the results 
of numerical simulations are not consistent. Several au- 
thors were led by their investigations to the conjecture 

' that Du = 4 ||7|,p^. On the other hand comparable sim- 
ulations in various dimension display no mean-field be- 
havior for D = A which was interpreted as an evidence 
that the values of the upper critical dimension is greater 

; than four 

In this paper we consider the BTW model in various 
dimensions and improve the accuracy of the analysis sig- 
nificantly. Our analysis reveals that the scaling behavior 
of the four dimensional model is characterized by the 
mean-field exponents with additional logarithmic correc- 
tions. We benefit in our analysis from the exact solution 
of the directed BTW model which displays logarithmic 
corrections at the upper critical dimension D„ = 3 p4[ . 
This solution is used in order to develop a scaling anal- 
ysis for the directed BTW model which takes these log- 
arithmic corrections into account. This type of scaling 
analysis will then be applied to the usual BTW model. 
The important result of this analysis is that the scaling 



behavior of the probability distributions as well as the 
usual finite-size scaling ansatz are affected by logarithmic 
corrections for D = A. These logarithmic corrections are 
a particular feature of the four-dimensional system and 
could not observed in higher dimensions. This suggest 
that the value of the upper critical dimension is indeed 
four. 



II. THE BTW MODEL 

We consider the ZJ-dimcnsional BTW model on a 
square lattice of linear size L in which integer variables 
-Er - represent local energies. One perturbes the system 
by adding particles at a randomly chosen site r according 
to 

Er 1-^ Er + 1. (1) 

A site is called unstable if the corresponding energy E^ 
exceeds a critical value Ec, i.e., if E^ > Ec, where Ec is 
given by Ec = 2D. An unstable site relaxes, its energy 
is decreased by Ec and the energy of the 2D next neigh- 
boring sites are increased by one unit, i.e., 

Er ^ Ey- — Ec (2) 

En-[i,r Enn,T + 1- (3) 

In this way the neighboring sites may be activated and 
an avalanche of relaxation events may take place. The 
sites are updated in parallel until all sites are stable. 
Starting with a lattice of randomly distributed energies 
Eg {0, 1, 2, Ec—l} the system is perturbed according 
to Eq. (|^) and Dhar's 'burning algorithm' is applied in or- 
der to check if the system has reached the critical steady 
state Usually one studies several different quantities 
in order to characterize the avalanches: the number of 
relaxation events s (size) , the number of distinct toppled 
lattice site a (area or volume), the duration t, and the 
radius r. In the critical steady state the corresponding 
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probability distributions should obey power-law behavior 
characterized by exponents r^, r^, Tj and according to 

p^{x) ~ (4) 

with X G {s,a,t,r}. Because a particular lattice site 
may topple several times the number of toppling events 
exceeds the number of distinct toppled lattice sites, i.e., 
s > a. It is known that multiple toppling events can be 
neglected for Z? > 3 i.e., the distributions Ps(s) and 

Pa{a) display the same scaling behavior and especially 

Scaling relations for the exponents Ts,Ta, Tt and can 
be obtained if one assumes that the size, area, duration 
and radius scale as a power of each other, for instance 

t - r'^*'-. (5) 

The transformation law of probability distributions 
Pt{t)dt = P,.(r)dr leads to the scaling relation 

Itr = — \- (6) 

The scaling exponents ^^x' are important for the descrip- 
tion of the avalanches properties and their propagation. 
For instance the exponent indicates if multiple top- 
pling events are relevant {'-^sa > 1) or irrelevant (7sa = !)■ 
Since the exponent 7^^ determine the scaling behavior of 
the avalanche area with its radius, ^ar is an appropriate 
tool to investigate whether the avalanche shape displays 
a fractal behavior or not. Finally, the exponent 74,- is 
usually identified with the dynamical exponent z. 

The measurement of the probability distributions and 
the corresponding exponents [Eq. (Q)] is affected by the 
finite system size L. If the avalanche exponents exhibit 
no system size dependence the finite-size scaling analy- 
sis could be applied In that case the probability 
distributions obey the scaling equation 

Px{x,L) = L-^- g^ixL-'^'^), (7) 

where the exponents have to fulfill the scaling equation 
llql. The exponent determines the cut-off 
behavior of the probability distribution and it was shown 
that Vx = Ixr (see for instance Q]). The advantage of 
the finite size scaling analysis is that it yields addition- 
ally to the avalanche exponents Tx the important scaling 
exponents: the avalanche dimension Va, the dynamical 
exponent Vt = z, etc. 

The value of the upper critical dimension £)„ of the 
undirected BTW model is not known rigorously. Several 
attempts were made to determine the value of Du using 
numerical simulations [[7P Jl^ , p^ . Usually one considers 
the probability distributions and compares the avalanche 
exponents with the known mean-field values (see for in- 
stance ]l6t). But due to the limited computer power 
the implementation of the higher dimensional systems re- 
duces considerably the system sizes L and consequently 
also the straight portion of the probability distributions. 
This makes a determination of the avalanche exponents 
via regression very difficult for _D > 3. This disadvantage 



can be avoided by applying a finite-size scaling analysis. 
Our results obtained in this way are consistent with the 
assumption that D„ = 4 and that the avalanche dimen- 
sion is = 4 for L) > 4 Q . 

Recently, Chessa et al. considered the BTW model in 
various dimensions using the same finite-size scaling anal- 
ysis Q. Compared to |^] they examined larger system 
sizes in _D > 3 and used an improved statistics (up to 10^ 
nonzero avalanches). From their results which differ for 
D > A from those in 0| they concluded that _D = 4 is not 
the upper critical dimension. Especially they obtained 
from their finite-size scaling analysis Va ~ 3.5 for the 
four-dimensional BTW model, i.e., the avalanches dis- 
play a fractal behavior already for D = 4. The origin of 
these conflicting results is that the used statistics (2 x 10^ 
nonzero avalanches) in is not sufficient. Especially 
the fiuctuating data points at the cut-off of the distribu- 
tion Pa{a) lead to uncertain results (see Fig. 5 in Q). For 
instance it is possible to obtain with this data a collapse 
of the distributions Pa{a,L) for values of the avalanche 
dimension between Va = 3.4 and Va =4.1. 

Thus, there is no agreement in the literature on the 
behavior of the BTW model in different dimensions: 
Chessa et al. concluded from their analysis that the up- 
per critical dimension is larger than four and that the 
avalanches displays fractility already for D = 4. On 
the other hand there exist several theoretical approaches 
which leads to the conclusion that = 4: Real space |^ 
as well as momentum space renormalization group 

analysis predicted both Du = 4. From their exact solu- 
tion of the BTW model on the Bethe lattice Majumdar 
and Dhar concluded that £>„ > 4 because the fractal di- 
mension of avalanche clusters must be lower than that of 
the embedding space ||l^. This leads the authors to the 
conjecture that the avalanches are compact for D < Du 
and fractal above the critical dimension. This fractility 
of the avalanche structure was already observed. Con- 
sidering the avalanche propagation in higher dimension 
it was found that the avalanches are characterized by a 
compact activation front for D = 3 and D = 4. For 
D > 4 the compact shape of the activation front is lost 
and several branches propagate through the system with- 
out coalesce together again (see Fig. 8 in Q]). Here, the 
avalanche propagation can be described as a branching 
process which is the main feature of the mean-field solu- 
tion of sandpile models (see for instance (l^). Assum- 
ing that the clusters are compact and neglecting mul- 
tiple toppling events (which is justified for D = 3 and 
D — 4 0,|l^) Zhang derived in the continuum limit the 
equation Ta — 2 — 2/D [|l8| which gives the mean-field 
value Ta —3/2 again for _D = 4. 

This incoherent picture of the behavior of the BTW 
model in higher dimensions leads us to reconsider the 
avalanche distributions again and compare our results 
with those of |^ . In contrast to our previous work Q we 
use now larger system sizes {L < 128 for D — 4, L < 48 
for D — 5, and L <24 for D = 6) and increase the statis- 
tics significantly, i.e., we averaged all measurements over 
at least 5 x 10^ non-zero avalanches. As usual we mea- 
sure the avalanche distributions [Eq. ^] by counting the 
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numbers of avalanches corresponding to a given area, du- 
ration, etc. and integrate these numbers over bins of 
increasing length (see for instance fl^). Successive bin 
length increases by a factor b > 1. Throughout this work 
we performed all measurements with the factor b — 1.2 
since larger values of b may change the cut-off shape of 
the distributions. Applying the finite-size scaling analysis 
this effect could lead to uncertain results for the scaling 
exponent Va (we found that this effect has to be taken 
into consideration at least for b> 1.5). 

We focus our attention on the finite-size scaling analy- 
sis. Performing this analysis it is informative to produce 
the data-collapse not only for all curves corresponding to 
different system sizes but also to check the obtained data 
collapse for selected curves. For instance the finite-size 
scaling analysis of two curves corresponding to two suc- 
cessive system size (Li < L2) allows to check whether 
the actual scaling regime is already reached. This anal- 
ysis is shown in Fig. |l| for the three-dimensional BTW 
model where it is known that finite-size scaling works for 
L > 64 . If one performs the finite-size scaling analysis 
for two system sizes with Li < L2 < 64 it is possible 
to obtain a data collapse (with small but systematic de- 
viations, especially at the cut-off) but then the scaling 
exponents depend on the system sizes. In Fig. |l| we plot 
the scaling exponent Va{Li, L2) as a function of the aver- 
age system size L = (LiL2)^^^- With increasing system 
sizes the exponent tends to the value Va = ?>■ For L>QA 
no significant system size dependence could be observed, 
i.e., a cross-over to the actual scaling regime where finite- 
size scaling works takes place at Leo ~ 64. 

Analogous to the three-dimensional model we perform 
the same analysis for D = A, D = 5, and D = 6 and 
plot the obtained results in Fig. |^. It is remarkable that 
within the error-bars the values i^aiD = 3)+l, I'aiD = 5), 




L 



FIG. 1. The finite-size scaling exponents i^a of the BTW 
model for various dimensions. The values of the exponents 
are obtained from the finite-size scaling analysis [Eq. (Q)] 
of two probability distributions corresponding to two differ- 
ent system size Li and L2 and are plotted as a function of 
L = {LiL2y^^ ■ In order to compare the different dimensions 
i^a + 1 is plotted for D = 3. 



and i^aiD = 6) display for small system sizes the same 
finite-size dependence whereas the behavior of the four- 
dimensional system differs significantly from the other 
dimensions. The conjecture that the system size depen- 
dence of Va is independent of the dimension (except of the 
case D = 4) implies that the cross-over to the actual scal- 
ing regime takes place at a comparable value Leo ~ 64. 
This could explain why the finite-size scaling analysis 
performed by Chessa et al. for _D > 5 yields exponents 
which are lower than the mean- field value Va = 4 
Their considered system sizes for D>5 are outside the 
scaling regime where finite-size scaling works. Thus, they 
(and of course all other previous numerical investiga- 
tions |0,|l^,^) observed only the cross-over to the real 
scaling regime and not the real scaling behavior itself. 

The significantly different behavior of the scaling expo- 
nent Va for £) = 4 (see Fig. ^ is remarkable since with in- 
creasing system size no cross-over to a scaling regime with 
a system size independent exponent Va could be observed. 
It seems that the scaling behavior of the four-dimensional 
model differs in principle from all other dimensions. A 
possible explanation is that the value of the upper critical 
dimension is £)„ ~ 4. Then the unique behavior of the 
exponent Va and the observed deviations to the expected 
pure mean-field scaling behavior for I? = 4 ||] could be 
explained by additional logarithmic corrections which af- 
fect the scaling behavior and which typically occur at the 
upper critical dimension. 

In the rest of this paper we will show that our re- 
sults are consistent with the assumption that the scaling 
behavior of the four-dimensional BTW model is char- 
acterized by the mean-field exponents with additional 
logarithmic corrections. In the next section we consider 
the BTW model with a preferred direction of the dy- 
namics. This directed BTW model is exactly solved 
and it is known that logarithmic corrections occur for 
= 3 0. The directed BTW model is therefore a 
suitable paradigm to learn how the logarithmic correc- 
tions enters the scaling behavior at the upper critical di- 
mension. This method of analysing will then be applied 
to the four dimensional BTW model in the last section. 



III. THE DIRECTED BTW MODEL AT THE 
UPPER CRITICAL DIMENSION 

In this section we consider the directed version of the 
BTW model which was introduced and exactly solved in 
all dimensions by Dhar and Ramaswamy iQ. Directed 
models are characterized by a preferred direction of the 
toppling rules. For instance a relaxation process takes 
place in a two-dimensional model if the energy of a given 
lattice site (i, j) exceeds the critical value Ec = D: 

— > Ei+ij + Ec/D (8) 
Eij+i Eij+i + Ec/D. 
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FIG. 2. Snapshots of three arbritrary chosen avalanche of 
the directed BTW model at the upper critical dimension for 
L = 128. 



One usually considers in simulations directed systems on 
a square lattice with periodic boundary conditions in the 
direction perpendicular to the preferred direction and 
open boundary conditions parallel to the preferred di- 
rection. The system is perturbed on the first line only 
(top of the pile) and particle could leave the system only 
on the last line (bottom of the pile) . 

Caused by the definition of the toppling rules no mul- 
tiple toppling events can occur (=> Ts = To). Since the 
perturbation takes place only on the top of the pile the 
average flux of particles through a surface in a given dis- 
tance from the top is constant. This flux conservation 
leads to the scaling relation 



1 

2 . 

Tt 



(9) 



According to Dhar and Ramaswamy the avalanche expo- 
nents of the two-dimensional model can be obtained by 
mapping the avalanche propagation onto a random walk 
and one gets Tq = 4/3 and Tt = 3/2, respectively. 

For D>i the exponents equal the mean- field values, 
i.e., Ta — 3/2 and Tt — 2, and additional logarithmic 
corrections to the power-law behavior occur in D = 3 
which is the value of the upper critical dimension . A 
snapshot of several avalanches of the three-dimensional 
model are shown in Fig. The shape of the avalanches 
reminds of a branching process which characterizes the 
avalanche propagation in the mean-field solution. 

According to the exact solution of Dhar and Ra- 
maswamy the mean square flux m{t) out of a given sur- 
face t is given by 



m{t) = Fit') 



(10) 
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FIG. 3. The probability distributions P{t > T) of an 
avalanche of duration greater than or equal to T for the di- 
rected BTW model at the upper critical dimension D„ = 3. 
According to Eq. (0) T P{t > T) is plotted as a function of 
InT. The dotted line is plotted to guide the eye. 



with F{t) - 1/lnt for £) = 3 IJ. Since the average 
flux through a surface t is constant in the steady state 
the probability distribution of an avalanche of duration 
greater than or equal to T scales in leading order as 



P{t > T) 



1 



m(T) 



InT 



(11) 



The corresponding plot of the rescaled distribution 
T P{t > T) as function of the duration in a logarithmic 
diagram is shown in Fig. ^ and confirms Eq. ([ll|). The 
scaling behavior of the corresponding density distribu- 
tion Pt{t) is then given by 



(12) 



In Fig. U we plot the rescaled distribution Pt (t) / In t as a 
function of the duration t. The rescaled distribution ex- 
hibits a power-law behavior with the exponent Tt = 2, in 
agreement with Eq. (|l^) . The inset of Fig. ^ shows that a 
fit of the unsealed distribution Pt (t) leads to lower values 
of the exponents Tt, i.e., a simple regression analysis can 
lead to the wrong result that the probability distributions 
are not characterized by the mean- field exponents. 

The scaling behavior of the average avalanche dura- 
tion (t) L confirms the relevance of the logarithmic cor- 
rections. Using Eq. (|l^) the average duration is given 
by 



{t)L = / tPt{t)dt ^ {\nLf 



(13) 



because in directed models the maximum value of the 
avalanche duration <max equals the system size L. The 
scaling behavior of the average duration clearly displays 
the relevance of the logarithmic corrections since without 
these corrections the average duration scales as ^ InL. 
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FIG. 4. The probability distributions Pt{t) of the directed 
BTW model for D = 3. According to Eq. Pt{t)/\nt 
is plotted as a function of the duration t. The dashed line 
corresponds to a power-law with the exponent rt = 2. In 
the inset we plot Ptit) vs. t. Here, the curvature is caused 
by the logarithmic corrections. The inset shows that wrong 
results of the exponents were obtained if one does not take 
the logarithmic corrections into account. 



In order to confirm this result we plot in Fig. |^ the square 
root of the average duration as a function of the system 
size L in a logarithmic diagram. The scaling behavior of 
the average duration agrees with Eq. (|l|). 

The scaling behavior of the probability distribu- 
tion Pa (a) of the avalanche area displays also logarithmic 
corrections. The area a of an avalanche of total duration t 
is determined by the average number of toppling events 
in each surface t' <t (see ||lj] ) and one gets to leading 
order 



a{t) = J2 "^(^) 



t'=i 



hit' 



(14) 



Instead of the usual scaling behavior a ^ which is valid 
for D > Du the leading order of the area scales with 
the duration as a ^ i^/ Int. Since the maximum value of 
the duration tmax equals the system size the maximum 
avalanche area scales as 



(15) 



The maximum area Omax determines the cut-off behav- 
ior of the probability distribution and Eq. ( |l5|) indicates 
that the usual finite-size scaling ansatz [Eq. (|^)] has to 
be modified in the presence of logarithmic corrections. 
In the following we derive this modified finite-size scal- 
ing ansatz which describes the scaling behavior of the 
avalanche distribution Pa{a) for D — D^- We assume 
that the leading order of the probability distribution of 
the avalanche area is given by 



Pa{a) 



a'" 



(16) 



FIG. 5. The square root of the average avalanche distribu- 



\l/2 



VS. L of the directed BTW model for D 



tion (t) ^ 

solid line corresponds to a logarithmic dependence of (t)^' 
according to Eq. (M) 



The 

1/2 



with the mean-field exponent = 3/2 and where the 
exponent of the logarithmic corrections Xa has to be de- 
termined. Comparing Eq. ( |l2| ) with Eq. ( |l6| ) the corre- 
sponding exponent of the duration distribution is given 
by xt = 1. Using the transformation law for probabil- 
ity distributions Pa{a)da = Pt{t)dt one can derive the 
exponent Xa- Inserting Eq. (14) into Eq. (nff) one gets 



PaHt)] ^ <-3(lni)-"+3/2 2- 



Inlni 
Int 



and analogous 



da(t) 
dt 



(lnt)2 ■ 



(17) 



(18) 



The term of leading order in the transformation law has 
to vanish and thus we get Xa = 1/2. 

Due to the logarithmic corrections of the probability 
distribution Pa (a) the simple finite-size scaling ansatz 
Eq. (|^) does not work. The simplest ansatz is to as- 
sume that the rescaled distribution Pa(a)(lna)^'' obeys 
the finite-size scaling equation 



Pa(a,L) (In a) 



(19) 



with the universal function g and where the scaling func- 
tion f{L) has to be determined. For low values of the ar- 
gument of the universal function (a <^ flmax) the rescaled 
probability distribution is independent of the system size 
and characterized by the power-law behavior g{x) x~'^°- 
only. Thus we obtain /(L) ~ «max- Using the known 
scaling behavior of Omax we get the modified finite-scaling 
ansatz 



Pa{a,L) (lna)-» 



L 



(lnL)^» .g(aL"2lnL). (20) 



We present the corresponding scaling plot in Fig. ^. The 
data collapse of the different curves corresponding to dif- 
ferent system size L confirms the above analysis. 
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FIG. 6. The modified finite-size scaling plot of the 
probability distribution Pa (a) of the directed BTW model 
for D = 3. The data collapse of the different curves corre- 
sponding to different system sizes L confirms Eq. (^o|). The 
dashed line corresponds to a power-law with the mean-field 
exponent Ta — 3/2. 

In summary we showed that the scaling behavior of 
the directed BTW model at the upper critical dimension 
is characterized by strong logarithmic corrections. These 
logarithmic corrections affect the usual probability dis- 
tributions [Eq. (H)], the scahng equations [Eq. (||)], and 
the finite-size scaling analysis [Eq. (^]. The corrections 
are relevant in the sense that one has to take them into 
account in order to describe the real scaling behavior, 
otherwise one gets wrong values for the exponents (see 
Fig. |). 

Additionally we simulated the directed BTW model 
for _D — 4 and performed a finite-size scaling analysis. 
In agreement with the exact solution of Dhar and Ra- 
maswamy the simple finite-size scaling ansatz, i.e., with- 
out logarithmic corrections, works quite well and the 
corresponding exponents equal the mean- field exponents. 
Thus, the logarithmic corrections to the scaling behavior 
occur only at the upper critical dimension Dc ~ 3 [M. 



(InL) 



and 



(InL)^o 



(22) 



The fifth introduced exponent describes how the area 
scales with the duration 



ait) 



(inty 



(23) 



At the upper critical dimension the avalanche and scaling 
exponents equal the mean-field values Tq = 3/2, rt = 2, 
i^a = 4, j/j = 2, and jat = 2. In this way the loga- 
rithmic corrections are determined by five non-negative 
exponents which have to fulfill two scaling relations. The 
transformation law of probability distributions leads to 
the first scaling equation: 



dt 



Tat 

~2~ 



(24) 



where we make use of the equation jat — ('''t — l)/(''a ^ 1) 
and assume that the term of leading order of the logarith- 
mic corrections has to vanish. Under the condition that 
the scaling behavior of the leading order of the maximum 
avalanche area amax is given by Eq. ( p2| ) and Eq. ( p3| ) we 
obtain the second scaling relation 



ttmax - a(imax) - (J^ 2.)^t7at+r„ 
Na = Tat + latNt. 



(25) 



Here, we use that the scaling exponents equal the 
avalanche dimension [va = lar) and the dynamical ex- 
ponent [vt = Jtr), respectively. The relation jtrjat = 
lar (l^ leads to Eq. (^5|). Thus, the logarithmic cor- 
rections to the usual scaling behavior are determined by 
only three independent exponents. 

Corresponding to the directed BTW model for D — £)„ 
we assume that the distributions obey the finite-size scal- 
ing ansatz 



Pt{t,L) {\nty 



ft{L) gt{t/t 

max J 1 



(26) 



IV. THE UNDIRECTED BTW MODEL FOR D = 4 

In the following we return to the investigation of the 
undirected BTW model for D — 4 and show that the 
avalanche distributions are characterized by logarithmic 
corrections comparable to the directed BTW model at 
the upper critical dimension. First we generalize the scal- 
ing equations by introducing certain exponents which de- 
scribe the logarithmic corrections. Guided by our previ- 
ous analysis we assume that the probability distributions 
are given by 



{inty 



and Paio) 



(Ina)^" 



(21) 



ft - a 

The maximal avalanche duration and area which deter- 
mine the cut-off behavior of the corresponding distribu- 
tions should scale with the system size as 



Pa{a,L){lna)-^- ^ fa{L)ga{ala^^^), (27) 
where the scaling functions ft and fa are given by 

ft{L) ^ i-^*-' (InL)^'^S (28) 



fa{L) 



L- 



(InL) 



(29) 



since for small values of the argument of the universal 
functions {t ^ tmax and a ^ a^^x) the probability distri- 
butions are independent of the system size. For D = D„ 
the avalanche and scaling exponents equal the mean-field 
values Ta = 3/2, rt = 2, Va = 4, and vt = 2. Thus the 
finite-size scaling analysis of the area and duration distri- 
bution gives the correction exponents Nt, Xt, Na, Xa, and 
the obtained values have to fulfill the equation [Eq. (12? 
and Eq. (I2F 
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FIG. 7. The system size dependence of the average 
avalanche duration (t)_L of the BTW model for D = A. To 
guide the eye we plot the solid line which corresponds to 
Eq. dsil). 



2{xt — Xa) 



2N,. 



(30) 



Analogous to the analysis of the directed BTW model 
we consider the scaling behavior of the average avalanche 
duration (i)^ before we apply the finite-size scaling anal- 
ysis. According to Eq. ( pi|) and Eq. (|2|) the average 
duration is given by 



{t)L = / tPt{t)dt 



In 



(Ini) 



Aft 



xt + 1 



(31) 



which allows additionally to the finite-size scaling analy- 
sis an independent determination of the exponents xt and 
Nt- We tried several values of Xt and Nt and obtained 
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FIG. 8. The modified finite-size scaling plot of the prob- 
ability distribution Pt{t) of the BTW model for D = 4 
and L = 24, 32, 40, 48, 56, 64, 72, 80, 96, 128. The data col- 
lapse of the different curves corresponding to different system 
sizes L confirms Eq. (p^). The dashed line corresponds to a 
power-law with the mean-field exponent Tt — 2. 



FIG. 9. The modified finite-size scaling plot of the prob- 
ability distribution Pa{a) of the BTW model for D = 4 
and L = 24,32,40,48,56,64,72,80,96,128. The data col- 
lapse of the different curves corresponding to different system 
sizes L confirms Eq. (^). The dashed line corresponds to a 
power-law with the mean-field exponent Ta — 3/2. 



a good result for « 1/2 and Nt 1/2. In Fig. ^ we 
present in an logarithmic diagram (t)^/'^ vs. L'^ /{\n L)^/^ . 
The plotted values are located on a straight line, in agree- 
ment with Eq. (|l|). 

These values are confirmed by the finite-size scaling 
analysis of the duration distribution Pt{t) according to 
Eq. (p^). A satisfying data collapse is obtained for xt = 
1/2 and Nt = 1/2 as one can see in Fig. || 

Now we consider the finite-size scaling analysis of the 
area duration Pa{a) . Since the correction exponents have 
to fulfill Eq. (30) it must be possible to produce the 
data collapse of Pa{a.) by varying one parameter only 
if we assume that the above determined values Xt = 1/2 
and Nt — 1/2 are correct. We eliminated Na in the 
scaling ansatz [Eq. ( ^ and Eq. (|2|)] and varied the 
exponent Xa- An almost perfect data collapse is ob- 
tained for Xa ~ 1/4. The corresponding scaling plot 
is shown in Fig. and confirms the accuracy of the 
above determined exponents xt and Nf. In contrast to 
the three dimensional model, where the simple finite-size 
scaling works for L > 64 Q , the finite-size behavior of 
the four-dimensional model are governed by the logarith- 
mic corrections and the modified finite-size scaling ansatz 
works very well already for L > 24. Finally we mention 
that Chessa et al. who used the simple finite-size scaling 
ansatz obtained a less accurate data collapse for L > 48 |§] 
since they did not take the logarithmic corrections into 
account. 



V. CONCLUSIONS 

We studied numerically the dynamical properties of 
the BTW model on a square lattice for Z? > 3. Our in- 
vestigation of the avalanche distribution which includes 
a careful examination of the finite-size corrections shows 
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that analyses @,|| of the BTW model for D > 4 are not 
conclusive. Our results are consistent with the assump- 
tion that the scaling behavior of the four-dimensional 
BTW model is characterized by the mean-field exponents 
with additional logarithmic corrections. We provide nu- 
merical tests for the theoretically predicted logarithmic 
correction terms for the directed BTW model at the up- 
per critical dimension £)„ = 3. We introduce a refined 
finite-size scaling analysis which takes these logarithmic 
corrections into account. These logarithmic corrections 
occur for D — 4 only, strongly suggesting that the value 
of the upper critical dimension is four. To proof this 
definitely in our opinion it is necessary to show that the 
distributions of the five dimensional model are character- 
ized by the pure mean-field values. Unfortunately, due 
to the limited computer power it is at present impossible 
to consider the actual scaling regime for D — 5. 
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